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I. INTRODUCTION 



One of the most prominent and important features of hadronic final states produced in particle collisions is the jet 
structure, i.e. the presence of a small number of collimated groups of particles recoiling against each other. Quantifying 
the structure of hadronic final states in e + e~ annihilation in terms of event shapes allows for a direct comparison of 
experimental data with predictions in QCD and thus for precise measurements of parameters and stringent test of 
the theory Q. Event shape variables avoid direct association of particles to individual jets and calculate instead a 
single number that accounts for, and classifies the event according to its jet topology. Generally, the observables are 
constructed such that a value close to zero corresponds to a dijet event topology where two jets of energetic particles 
are produced back-to-back with additional soft particles between the jets. Most events enter in this dijet region of 
the event shape distributions because two jets can already be produced at tree-level while three and more jets are 
suppressed by additional powers of the strong coupling. 

Among the most common event shape variables are the thrust T defined by 

T ee max, (^0] , (1) 



Ei \Pi 

where the sum is over all particles, and pi is the three- momentum of particle i. The thrust axis is the unit vector n 
which maximizes the expression in the parentheses. Since T ~ 1 characterizes dijet-events it is convenient to use the 
variable 

T ee 1 - T (2) 

as the thrust variable. The plane through the interaction point and perpendicular to the thrust axis divides the event 
into two hemispheres, Hi and Hi. The maximum of the squared invariant masses M 2 and Mf of all the particles in 
Hi and Hi, respectively, defines the heavy jet mass variable || 

max(M 1 2 ,M|) 

P= qi , (3) 

where Q is the e + e~ cm. energy. In the dijet limit where r <C 1, the thrust can be written as a function of the 
hemisphere Masses M\ and Mi, 

In fact, the double-differential M\-Mi invariant mass distribution represents by itself an event shape distribution, 
where small values of M 2 2 /Q 2 correspond to the dijet region. 

For massless quarks in the dijet region the M\-Mi distribution can be described by a factorization theorem B, la, 

|, 0, 1 § 



a H Q (Q,n) / dl + d£-Ji{Ml - Q £ + , fx) J 2 (Ml ~QT,y) S '{£+,£-, ft), (5) 



^dijet 

dMfdMj 

which is valid at leading order in M 2 2 /Q 2 . Here, cto is the tree level total cross section, Hq is a calculable hard 
coefficient and J\ % i are calculable jet functions, whereas S(£ + ,£~,fi) is the hemisphere soft function. A similar 
factorization theorem for the M\-Mi distribution can be derived for the production of quarks with mass m 3> Aqcd- 
Here the dijet region corresponds to the region near the heavy quark mass resonance, where §i t 2 = (M 2 2 — m 2 )/m -C m. 
The corresponding factorization theorem has the form j|, ||] 

d 2 cr d » ot TT , TT ( Q \ /"„.,._„ /. Q£ + \ „ /. Q£ 



dM 2 dM 2 



°o Hq(Q, fi m ) H rn (m, —,(A m ,fA [d£ + d£ B + (s 1 - ^—,fx) B_ (s 2 - ^—,fx)S(£ + ,£ , (6) 
V m J J \ m / V ml 



which is valid at leading order in m 2 /Q 2 and (§i t 2 + Y q )/m. Here H m is a calculable hard coefficient and B± arc 
calculable jet functions for heavy quarks describing invariant mass fluctuations below the scale m. The term T q is 
the heavy quark width. The soft function S is equivalent to the one in Eq. (|^). As a consequence of the separation 
of the different physical modes in the factorization theorems (||) and (0), the hard coefficients, and the jet and soft 
functions are renormalization scale dependent. 
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The soft function carries information on how the soft radiation between the two energetic jets is associated to the 
invariant masses Mi and M2. For the hemisphere prescription described above it is defined as (J, |[ |[ 0, [l0| 

s(£+,r,ri = -^^^+-fc+ a )5(r-fc s - b )(o|(F fi ) cd (y n ) ce (o)l^}(x s |(y„t)^ (Fjo*(o)|o>. (?) 

Here kf a is the total plus- momentum of soft hadrons in X s that are in hemisphere 1, k~ b is the total minus momentum 
for soft hadrons in the other hemisphere. The soft function for thrust is related to the hemisphere soft function by 

St(t,/x) = J d£ + drs(r- e+ + £ }s(£+,r,fx). (8) 

The definition of the soft function only depends on the light-like kinematics and the color state of the primary 
quark-antiquark pair and thus can be written in terms of the Wilson lines 

Y^(x) = P exp {^9 j dsn-A s (ns+x)j , Y n \x) = P exp {^9 J dsn- A s {ns+x)^j . (9) 

In general, S(£ + ,£~,fi) is a nonpcrturbative function that peaks for £ ~ Aqcd when \i > Aqcd- Depending on 
the size of M^a different aspects of the soft function arc important since the convolutions in Eqs. (@) and (^) probe 
momenta ^ ± ~ M^ 2 /Q and ~ si^ra/Q, respectively. In the immediate resonance region we have Mf 2 ~ QAqcd and 
si,2 ~ QAqcd/™ + P<j H) 1 , and the nonperturbativc distribution described by S(£ + ,£~,a) affects directly the shape 
of the differential cross section. Here, the soft function can be written as a convolution of the partonic soft function, 
computed in perturbation theory at a scale jj, = ft Aqcdj with a nonpcrturbative model function that can be 
determined from experimental data [jfl] , 

p + 00 r+oc 

S(£ + ,£-,fi)= d£'+ d£>- S part (£+ -£'+,£--£'-, fi)S niod (£' + ,£'-). (10) 



In the tail region away from the resonance we have Mf 2 3* QAqcd and s\ t 2 ^ Q Aqcd/w + r ? and the soft function 
can be determined from the partonic soft function plus power corrections that can be also related to Eq. (|l^). Note 
that the partonic soft function contains S functions and plus-distributions of the variables £^ , so that the lower limit 
of the ^-integrations in Eq. ( fio| ) is zero. In Ref. the function 

was suggested as a two-parameter model for Sraod- Here Af(a,b) is a factor that is chosen such that the integral of 
/exp of the positive ^ /=t plane is unity. This normalization is required by the consistency of Eq. (^|) with the power 
expansion in the tail region. The 0(a s ) corrections to the partonic hemisphere soft function were computed in Ref. j^] 
(see also Ref. [|l2| for the 0{a s ) corrections to the partonic soft function for the thrust distribution). 

It has been noticed in Refs. that the partonic threshold of the soft function S^art (^ 7 m) at £ ± = has 

an O(Aqcrj) renormalon ambiguity that is similar in nature to the well-known O(Aqcd) pole mass renormalon. In 
Ref. [ [TT| is was shown that this ambiguity can be removed by introducing a gap in the soft function model such that 
it vanishes for £ /=t < A. A way to achieve this is by defining 

Smod^'V") = U P (£' + - A,f- - A) . (12) 

Here A can be interpreted as the minimum hadronic energy deposit in each hemisphere. Via the convolution in 
Eq. (0) the term A compensates the renormalon ambiguity in at its partonic threshold. It is then possible 

to explicitly remove this renormalon by writing A = A + 5A, where A is free of an O(Aqcd) renormalon and SA 
is a perturbation series that cancels the renormalon ambiguity in the partonic soft function order-by-order. It was 
demonstrated in Ref. [0 that this procedure leads to a substantial reduction of the perturbative uncertainties that 
come from the soft function. 



1 The heavy quark width Tq plays a vital role for top quarks, but can be neglected for bottom quarks. 
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In this paper we determine the O{ot 2 s ) two-loop corrections to the partonic hemisphere soft function S'part, which 
previously was the only term in the factorization theorems (0) and (^J) that was not known at 0{a 2 s ). The result enables 
a full next-to-next-to-leading logarithmic (NNLL) determination of the distributions described by the factorization 
theorems (|5|) and (^J). Using constraints on the form of the soft function from its rcnormalization group (RG) 
properties, the nonabelian exponentiation theorem |14, [l5| and from the 1 <-> 2 symmetry with respect to the 
hemispheres Hi and H2, it is possible to determine the analytic form of the previously unknown 0(a 2 s ) contributions 
of the soft function up to two constants. These constants can be determined from Eq. (|^) by taking the known 
two-loop results for the hard coefficient H q ]l6|, 17|, [l8| and the jet function Ji i2 |l9| as an input and using numerical 
results from the MC program EVENT2 ]20|, |21] for thrust and heavy jet mass distributions for massless quarks at 
0{a 2 s ). The result we obtain are confirmed by correct predictions for distributions of variants of the thrust and heavy 
jet mass variables that can also be obtained from EVENT2. To our knowledge these variants of the thrust and heavy 
jet mass variables have not been defined in the literature before. 

The program of this paper is as follows. In Sec. || we review the previously known properties of the partonic 
soft function. We show that these constraints determine the two-loop partonic soft function up to two unknown 
constants that can be determined numerically. In Sec. |l| we determine the 0{a 2 ) expressions for the different event 
shape distributions employed in our numerical work and in Sec. |y| we present the analysis to determine the two 
constants from EVENT2. In Sec. |v| we introduce a renormalon-free gap parameter A from the position space soft 
function. This gap parameter depends on an infrared scale R and the renormalization scale fx, and we determine its 
evolution equations in R and /i. Section VI contains a brief numerical analysis illustrating the impact of using the 
renormalon-free gap parameter and of accounting properly for its scale dependences. The conclusions are presented in 
Sec. VII. We have attached three appendices collecting useful formulae on the Fourier transform of plus-distributions, 
on results for the hard coefficient Hq and the massless jet function J adapted to our notation, and on the cumulative 
event shape distributions that are used in this work. In particular, App. [A| contains a determination of the previously 
unknown non-cusp NNLL anomalous dimension of the hard coefficient H m that appears in the factorization theorem 
for massive jets in Eq (H). 



II. PROPERTIES OF THE HEMISPHERE SOFT FUNCTION 

In this section we summarize the known properties of the hemisphere soft function. These properties lead us to a 
particular analytic form for the previously unknown parts of the partonic soft function at 0{a 2 ) that depend only 
on two unknown parameters. These parameters are determined numerically in Sec. IV. To simplify the notation 



and avoid cluttering due to convolution integrals in the variables we use in this work mainly the position space 
representation. A number of formulas in the £^ momentum space variables useful for future applications are collected 
in the appendix. In position xi^-space the soft hemisphere function is defined as, 



dTdt 



S(t + ,£-,ti). 



(13) 



The soft function is x\ «-> X2 symmetric because the definition of the soft function in Eq. (0) is symmetric with 
respect to exchanging the hemispheres H\ and H 2 . 

Result at 0(a s ). It is straightforward to compute the partonic hemisphere soft function at 0(a s ) from the definition 
given in Eqs. (|7|). The result reads M 



Spaxt {Xl,X 2 ,fl) = 1 



In (i xie' <E fi) + — 
o 



C F a s (/j,) 



In 2 (ix2e lE n) + — 



(14) 



Renormalization Group Structure. From the dijet factorization theorem for massless jets in Eq. (^]) one can derive 
consistency conditions |^ which relate the RG-evolution of the soft function to the RG-evolution of the hard coefficient 
Hq and the jet function J. Details of the computation can be found in App. |aJ. Given S(xi , X2, £*o) at the scale \i = (1q 
we find that S 



S(xi,X 2 ,H) = U s {x 1 ,^,fi )U s (x 2 ,^^o)S{xi,X2,^o), 



where 



U s (x, (a, no) = exp 



;(r s ,/i,(U ) In (ixnoe lE ^j + A'(L S , 7 S , /i, ^ ) 



(15) 
(16) 



5 



with [r = a 8 (^i)/a s (no)] 



Q(T a ,fJ,,fJLo) = - 

Po 

~ (Is 



ft\ a s (/i ) 



4tt 



(r-1) + 



K(r s ,is, n, no) =w — ,/x,Mo) - 



47T 



4tt 



r o 

ftr* ft , /ftr s ^ 



1 In r 

r 



rf _ ft _ A. 
r° ft ft, Vrg 

rj_ft 
r2 /3o 



'o 



(1 — r + lnr 



ftrg #, VA>r° ft 2 



(r — 1)1 lnr — 



PiTl ft 

ftr° ft 



2(4tt) 2 

+ A 
2ft 

(r-1) 



r 2 -l 



In 2 r 



ftr* ft r 
ftr° + ft r° pi) 2 



(17) 



where for uj the N fe LL solutions correspond to the terms up to order a*, and for K the N fe LL solutions correspond to 
the terms up to order o^ -1 . The solutions depend on the cusp and the non-cusp anomalous dimensions of the soft 
function and the beta- function [a s = a s (/x)], 



r s [ a s 
da* 



E 

fc=0 



4tt 



k+l 



r* 



J s [a s 



E 

fe=0 



a,\ k+1 



47T/ 



7 fe 
is > 



dln/i 



ft. 



fc=0 



ft, 



(18) 



The first few coefficients sufficient for NNLL order running read 

- r° = -4Cf 

X CUSp w ^ 5 



r° = 
ri = 

r 2 = 



cusp 



490 536, 44 4 88 

~ 1 cusp =^^1 — + -^7T - — 7T - — C 3 ) + U A L F ±n 



/1672 160 



V 27 



27 



224 
— C3 



220 



C F Tn f [ — ~M^)+-C F {Tn f ) 



64 



0: 

■ Ca Cf 



808 
~2T 



11 



28 Cj 



and we also have 
ft = 
ft = 



— Ca Tnr 

3 3 f 



11 

2857 
~54~ 



ft 



224 4 
"27~ ~ 9 



(19) 



34 



20 



C A Tn f - AC F Tn f . 



1415 2 158 x2_ L 44 r1 N2 

— — C A Tn/ + — C A (Tn f ) + —C F (Tn f ) 



205 
~9~ 



C A C F Tn } +2C%Tn f 



(20) 



for light flavors, from the running of the strong coupling in the MS scheme. 

It is an important feature of the factorization theorem that the RG properties with respect to the variables x\ and 
x 2 factorize such that in S(x\, x 2 , (J-) the //-dependence can only arise in terms of powers of ln(xi t2 fi). As shown in 
Ref. [pTj it is therefore possible to write Eq. (p~5h in the form 



S(xi,X2,n) = U s (xi,fi, {ixie< E ) J ) U s (x 2 , M, (i x 2 e lE ) x ) S{x\, x 2 ) , 
where S is /x-independent and S(xi,x 2 ) = S(x 2 , x\). 



(21) 



Nonabelian Exponentiation. In Refs. 14, it has been proven to all orders in perturbation theory that QCD matrix 
elements with arbitrary number of external gluons exponentiate, if their operator definition can be written entirely in 
terms of Wilson lines. If the external gluon final state is symmetrized, the exponentiation also holds for contributions 



6 



of these matrix elements to cross sections and production rates. The exponentiation property also applies to the 
hemisphere soft function and is most transparent in position space. Using Eq. ( pl| ) we can therefore write, to all 
orders in perturbation theory, 



?part 



( Xl ,x 2 ,fx) =U s (x 1 , l i,(ix 1 e^)- 1 ) U s (x 2 , f i,(ix 2 e^)- 1 ) e T ^' x ^ 



exp 



K(T s , ls ,fi,(ix 1 e /B )- 1 ) + K (T s , ls ,[i,(ix 2 e~ fE )- 1 ) +T(x u x 2 ) 



(22) 



where T(xi,x 2 ) — T(x 2 ,xi). In this context exponentiation means that the argument of the exponential has a 
simpler color structure than <S P art(£i, x 2 , /i) itself. There are two specific features that are worth to mention: (i) At 
0(a") the highest power of logarithms of x\ j2 in the exponent is ln n+ , while in 5* it is In 2 " and (ii) the exponent does 
not contain any a™C F terms except for n = 1. The latter property means that in QED the exponent is 0(a s )-exact 
and does not contain any higher order corrections. 

Analytic form of T{x\ 1 x 2 ). It is now straightforward to determine the analytic form of the function T{x\ 1 x 2 ) to 
0(a 2 ). If there is only a single gluon in the final state it is either in hemisphere 1 or in hemisphere 2. Thus to 0(a s ) 
the soft function is a symmetric sum of two terms each of which is either a function of x\ or of x 2 . If there are two 
partons in the final state, they can be in different hemispheres, and a non-trivial dependence on x\ and x 2 can arise. 
Thus to 0(a 2 ) the function T must have the form [/i Xl 2 = (ixie 7s ) -1 ] 



T(xi,x 2 ) = — - h H h + 2 j— g t 2 {xi,x 2 ) 

An An {Any 



(23) 



where t\ that can be read off from Eq. (14) 



(24) 



Moreover we have t 2 (x\,x 2 ) = t 2 (xi/x 2 ) = t 2 (x 2 /xi), since the x±_ 2 arc variables which have the dimension of an 
inverse mass. Here it is worth to mention that the term a 2 , actually reads c (a 2 (/i Xl ) + o. 2 (^i X2 )) + da s (/j, Xl ) a s {jjL X2 ) 
with 2 c + d = 1. Since any event shape distribution in the dijet limit can be written in terms of delta functions and 
plus distributions, we can further use the information that the momentum space soft function only depends on 
through 5-functions or plus-distributions. Thus in position space t 2 can only depend on x\ >2 through even powers of 
\n(xi/x 2 ), see Eqs. (C4. C5|). Using also the constraint that t 2 cannot contain any term ln"a;i,2 with n > 3, it must 
have the simple form 



t 2 {xi,x 2 ) 



si 



s 2 In 2 ( — 

\x 2 



(25) 



An important consequence of the exponentiation property is that si and s 2 do not contain any contribution 
proportional to the color factor C\. This feature serves as an important cross check for the numerical analysis we 
carry out in Sec. [V. 



Form of the Soft Function. Expanded in ce s (/i) and using the results collected above, the 0(a 2 ) position space 
hemisphere soft function reads [a s = a 8 {n),Li 2 = hi{i x\ y2 e lE ^)] 



^part 



(x 1 ,x 2 ,ijl) = 1 + 



Utt 



Ci 



( C F a t 



\ An 



-2 1 



ALi 



ALi 



1 + ~~2 



4L 2 + i 1 



8C 2 F {Li + L$) + 



- 2 -f + \A + ^C F Tn f 



— qCaCf 



(Ll + L 2 



32 



C F Tn 



(Ll+Ll) 



808 22 9 \ / 224 8 9 

" ~27 ~ ~9 n + 28Cs ) + CpTnf ( ~Tf + 9 



(L!+L 2 ) 



+ —C F + 2t 2 (x u x 2 ) } . 



(26) 
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The corresponding momentum space hemisphere soft function has the form [C± = l/^[0(£ ± ) ln™(£ ± //i)/(£ ± //x)] ] 



*part 



(t+j-, ^) = s(e + )S(r) 



Cp a s 
An 



8C\ + -6(£ + ) 



^C A C F - ^C F Tn f 



-8£l + ^6(i-) 



32C| 



536 



160 



6(e + )c L _ + 5{r)c\ 



- l2n z C F + C A C F ^- — + -tt z j + —C F Tn f 

( 808 22 9 \ 

64C 3 C| + C A C F f — - —n 2 - 28C 3 J + C F Tnf [ - — + -tt' 



40 

f 2i 2 (£+,r) 



134 



176 



- 77^ ^f + CaCM - «t 4 + — Cs + C F Tnt I - — - — 



27 



9 



9 



224 

"2T ' 9' 
40 



S(£+)C° 



64 



27 



9 



2 5(nw) 



(27) 



where 



t 2 (^+,0 = Sl 5{£ + )5{£-) + s 2 1 24 - y^ + ) 



<5(r 



2zCL - -<5(r 



-2£° £°„ 



(28) 



Note that while the distributions £J that appear in Eq. ( p8| ) depend on the renormalization scale /z, this /x-dependence 
cancels in the combination of all terms. For completeness we also present the RG evolution of the momentum space 
soft function. It has the form S(£ + ,£-,fi) = J d£'+£'-U s (£+ - £'+ , fi, IM>)U S (£- - £'~ , fi, fx )S{£' + ,£'~ , fx ) with 



U S {£,H,IM>) 



Mor(-w(r s ,/j,/i )) 



£i+di(r a ,/i,/io) 



(29) 



where Q and if are given in Eqs. (|T^). The result is obtained without any effort using the position space result in 
Eq. ( |l6| ) and the Fourier transformation given in Eq. ( ]C4| ). The plus function with an arbitrary exponent l + Q with 
uj < 1 is defined by 



(x) 



1+Cj 



lim 

0^0 



(30) 



III. THRUST AND HEAVY JET MASS DISTRIBUTIONS 

In this section we determine the 0{a 2 s ) fixed-order cumulative distributions for the thrust and the heavy jet mass 
variables, and for their generalizations, called r Q and p a and defined below. The distributions are derived from the 
dijet factorization theorem (^|) for the double differential hemisphere invariant mass distribution for massless quark 
production. At 0{a 2 s ) the distributions depend on the parameters s\ and s 2 in the hemisphere soft function that 
arc determined numerically in Sec. |v|. Since EVENT2 produces 0(a 2 s ) distributions for the renormalization scale 
fx = Q we also use this scale choice for all the evaluations that follow in this section. This means in particular 
that any summation of logarithms in the hard coefficient, the jet function or soft function is neglected. Since the 
renormalization scale is fixed to /i = Q it is convenient to use the dimensionless variables iri )2 and 

A>2 = ^fh ~ t ± = 
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We start by writing down the double differential hemisphere invariant mass spectrum in the dijet region shown in 
Eq. (H) in position space representation normalized to the tree level total cross section, 



1 da dijet 



= / dM{ dM£ e 



1 da d » et 



cr &M\ dMi 



do dxi dx 2 

= Hq(Q, Q) Ji(ii/Q 2 , Q) J 2 (x 2 /Q 2 , Q) S(xi/Q, x 2 /Q, Q) 
We now define the event-shape variables t q and p a as 

2 aM?+M$ _ 2 (a Mf + M%) 



(32) 



1 + a 
2 

T+a Q 2 



— max(a Mf, M|; 



1 



Q 



max(aM 2 ,M 2 2 ) 



(33) 



For a = 1, T a and p a reduce to the common thrust and heavy jet mass variables. For a^l the two hemispheres get 
different weights and probe the M 2 -M| distribution asymmetrically. From the M 2 <-> Mf symmetry one can derive 
the relations 



T l/a 



and 



Pl/a 



(34) 



We note that, if a is chosen much larger or smaller than one, the r Q and p a distributions in fixed-order perturbation 
theory develop large logarithms of a. These logarithms arc examples of non-global logarithms |2^], and they can be 
summed in the factorization theorem by an independent setting of the renormalization scales that govern the invariant 
masses of the two hemispheres. It is now straightforward to determine the 0(a 2 ) fixed-order cumulative r Q and p a 
distributions in the dijet limit for p = Q 



£^ jct (ft) = f 
Jo 



n , 1 dcr dijet 

dr a 

Co ar a 



dr a dM{ dMi d(r a 



o 



l + Q 



(aMi + Mi))- 



1 da A ^ et 



cr dM 2 dM 2 



2tt x — iO 



(35) 



rfl i j dijet /-H i 

Spl jet (r!) = / dp a = / dp Q dM 2 rfM 2 2 5(p Q - T ^max(aM 2 ,M 2 )) — 

JO CTO dp a J (T 



1 da diiet 



(T dp a 

+ °°dx\ dx 2 
_^ ~2tt ~2tT 



(to dMfdM^ 



(x\ — iO)(x 2 — iO) 



(36) 



The analytic expressions for the cumulative distributions above are given in App. |B| Note that the function t 2 in 
Eqs. (^6j) and ( p7| ) docs not lead to any logarithmic terms involving the cut-off Q in the cumulative distributions. 
While in one dimension such a behavior can, in momentum space, only be obtained from a delta function, in several 
dimensions, it can also be achieved by proper combinations of plus-distributions as shown in Eq. (|28|). This behavior 
is easier to see in position space where t 2 only depends on the ratio x\/x 2 and contributes as 



in Eqs. (|3 
value of a 



and I 



, 27T 

+ °°dx\ dx 2 



x — iO 

i fi (xi+x 2 ) 



si + s 2 In a 
t 2 (ax 1 /x 2 ) 



6(0.) 



si 



si In a 

3 



(37) 



, 2?r 2tt (xi - iO)(x 2 - iO) 

The constants depend on the definition of the event-shape variable and in particular on the 
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IV. NUMERICAL ANALYSIS USING EVENT2 



In this section we determine the constants s\ and s 2 in Eqs. (|25|), (|2q), which cannot be determined from the 
general arguments discussed in Sec. |l[ 

Method. The EVENT2 program determines numerical estimates for event-shape distributions at 0(a 2 s ) in the fixed- 
order expansion for fj, = Q in full QCD. Using the variable y gencrically for r, r Q or p ai the distributions have the 
form 



1 da 
cr dy 



My) 



a s (Q) 



2n 



The corresponding cumulative distributions read 



E„(n) 



1 da 



a dy 

Defining bin boundaries O™ for n = 0, 1, 



2tt 



q.(Q) 

2tt 



E«(f2) 



C(y) 



a.(Q) 
2tt 



with < 0" < f2" +1 and f2 T 



(38) 

4 2 ' (39) 
1, EVENT2 can determine the 



sum of weights of events falling into the n max bins, which for the nth bin represents a numerical estimate for 



Aa n = E v (n n ) - lyfr- 1 ) = / d y — ^- 

J nn -i a dy 



Aa^ 



2ir 



<*s(Q) 

2ir 



(40) 



To determine the unknown constants s± and S2 in the dijet distribution from Eq. (^]), one can use the fact that they 
do not appear in Aa n for any choice of bin boundaries. This is because a dependence on s li2 can only appear in 
integrations that contain the threshold at y = 0, sec Eqs. (|35|), ( |36| ) and (37). The method starts by subtracting the 
known dijet contributions Acr„ ' dl ' et from the full theory Aa^ computed by EVENT2. For da/dy the difference is 
at most logarithmically singular for y — > and thus intcgrablc at y = 0. Thus for the cumulative distribution the 
remainder 



E^ cst (r>) = e„(0) - Ef ot (r!) 



(41) 



vanishes for CI — > 0. The remainder distribution is also independent of si and s 2 for any fi > 0. It is the aim 
to determine a numerical estimate for £( 2 '' rest (l) from EVENT2. Using that T, y (l) is equal to the total cross 



section 



E„(l) - R(Q 2 ) 



Otot 



C F a s (Q) 
2ir 



a s (Q) 
2ir 



i 2 



Cp + Ca Cf 



123 




11 1 / 


11 Ca 

8 s 


+ C F Tn f 


-y + 4 C3 _ 



(42) 



we can then determine the constants ti.y for the different event-shapes using Eqs. 

E^' dijot (l) = ^ _ sCaJ.reet^j _ 



) and (B3) and the relation 

(43) 



Sum of all Color Factors. For our numerical analysis we ran EVENT2 with 10 10 events for rif = 4 and for rif = 5. 2 
We used 80 bins from 10 -4 to 1 with logarithmic bin boundaries located at Q n = lO'" -80 )/ 20 with n = 0, 1, . . . , 80. 
Since EVENT2 works with an internal infrared cut-off (set to 10 -8 ), it is not possible to obtain numerical results for 
a bin with a lower boundary located at fl = 0. Thus to obtain a numerical estimate for Ey ost (l) one has to rely on 
an extrapolation of the numerical results for Ey 2 '' rGSt (l) — Ey 2 '' rcst (f2) taking the limit fi — ► 0. 



2 The task was distributed over 100 parallel jobs and took less than 30 hours to complete for each value of rif. 
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FIG. 1: 



Results of the numerical analysis to obtain Ey 2 ' ) ' rcst (1) and t2, y for y 



t (top), n (middle) and pi (lower row) for n/ = 5. 
The first column shows the full QCD results for the ^-distributions in bins with center at Q obtained from EVENT2. The 
second column shows the QCD distributions minus the singular dijet contributions described by the factorization theorem (||) . 



The third column shows E 



(2), rest 



(!)-£ 



(2), rest 



(SI), and the last column shows tz,y, taking E 



(2), rest 
V 



(1)-E 



(2), rest 



(n) as an estimate 



for E^ re8t (l). The solid and dashed lines in the panels in the column on the right represent the central value and the errors 
for t2, y which are estimated from the limit of small S7. 



Figure [T] shows the results of our numerical analysis to obtain estimates for T,y s (1) for the variables y = r, t\ 
and pi for rif = 5. 3 The first row is for t, the middle row for n and bottom row for p\. We note that the thrust r 
and the variable T\ agree in the dijet limit, but they differ concerning power corrections and in the multi-jet region 

(2) 

where r ~ n ~ 0(1). The panels in the first column show the binned full QCD distributions Aa^n as obtained from 
EVENT2. The central value for each bin is displayed as a (colored) symbol and the statistical error as a vertical line. 

( 2 ) rest 

The panels in the second column show the remainder distribution Acr^n after the singular dijet contributions have 
been subtracted. We see that the remainder distribution falls off to zero for small SI and that the asymptotic regime 
O — > appears to be reached already for SI < 10~ 2 ' 5 . The statistical errors grow for decreasing SI for the remainder 
distribution because EVENT2 attempts to obtain a constant relative statistical error for a given binning in the full 



3 Our results for rij = 4 are analogous and not discussed in detail. 
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FIG. 2: Tree-level (green dotted line), one-loop (red dashed line) and two-loop (blue solid line) soft function S(£,£,fi) (see 
Eq. (p"o|)) without renormalon subtraction for fj, = 1.5 GeV (a s (1.5 GeV) = 0.3285) and nj = 5. For the model function the 
parameters A = 0.55 GeV, A = 0.1 GeV and (a, b) — (3.0, —0.5) are used. 



QCD distributions shown in the first column, and because the contributions from the singular dijet terms increasingly 
dominate in size for £7 — > 0. The panels in the third column show the integral over the remainder distribution 



-,(2),rest 



(1) 



-,(2),rest 



display the estimates for the constants t^ y adopting £y rest (l) 

-,(2),rest 



,(2),rest 



(fi), which has to be extrapolated for O — » 0. Finally, the panels in the column on the right 

(Q) as shown in respective panels in the 

third column for £jf MeB ''(l) in Eqs. (|43|) . We use the average of the central values for log 10 fl < —2.5 to estimate 
our final numbers for t2, y and adopt the error at log 10 Q w —3.25 as the uncertainty. The results of the estimate are 
illustrated by the solid and dashed horizontal lines in the panels in the last column. We note that our results arc fully 
compatible with the theoretical expectation that t% T = t2,n- Using the relations 



si 



1 



{t2,T + *2,Tl) 



S-2 



{t-2,T + ^2,ri) — t 



2, Pi 



(44) 



we obtain 



si = 



-39.1 ±2.5 
-53.3 ±2.5 



5) 
4) 



S2 = 



-15.4± 1.5 
-14.9 ± 1.5 



5) 
4) 



(45) 



A method similar to ours has been applied recently in Ref. |26| to determine 0{a 2 s ) corrections to the soft function 
for thrust given in Eq. (|§|). The thrust soft function depends on si, but has no dependence on S2, see Eq. (|37j). 
Transferred to our notation their results reads S\(nf = 5) = —40.1 ± 3.1 and S\{rif = 4) = — 54.4 ± 3.0. The results 
are compatible to ours. In Ref. |2(| 10 10 events were used, but with linear binning and for > 10~ 3 . 

In Fig. U the tree-level (green dotted line), 0(a s ) (red dashed line) and 0{a 2 s ) (blue solid line) soft function 
S{t + ,t~,ijL) including the convolution with the model function as defined in Eq. ( [To| ) but without renormalon sub- 
traction are plotted over t = for fj, = 1.5 GeV (a s (1.5 GeV) = 0.3285) and .rtj = 5, and using A = 0.55 GeV, 
A = 0.1 GeV and (a, b) = (3.0,-0.5) for the model function in Eqs. (|ll|) and @. For the partonic 0{a 2 s ) soft 
function we have adopted the results in Eqs. J45| ) for the constants si,2. The corresponding uncertainty in the 0(a 2 s ) 
soft function is visualized by the additional light blue solid lines. The uncertainty is at the percent level where the 
soft function is large and absolutely negligible in comparison to the remaining perturbativc QCD uncertainties.. We 
will therefore adopt the central values given in Eqs. (|45| ) from now on without further discussion. The rather poor 
perturbative behavior of the curves shown in Fig. |^ with the unphysical negative values and the significant changes in 
the shape at higher orders is symptomatic for any choice of model parameters and renormalization scale and illustrates 
necessity to introduce the renormalon-free gap parameter. This will be discussed in Sec. 0. 



Cross Check using r a and p a for a ^ 1. An important cross check of the results in Eqs. (|45|) and also of the form 
for the function ti in Eqs. ( p5[ ) and ( p8| ) is provided by comparing predictions for ti, Ta and i2 iPa for a^l with the 
corresponding results obtained from EVENT2 based on a numerical analysis analogous to the one in the previous 
section. 
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FIG. 3: Analysis and results with errors for t2, y for y — r a (blue symbols and lines) and y = p a (red symbols and lines) with 
a = 1,2,3,5, 10 and 1/2, 1/3, 1/5, 1/10 and using n S = 5. 
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FIG. 4: Estimates for t2, y with uncertainties for y = r a (lower red dots with error bars) and y = p a (upper blue dots with 
error bars) obtained from EVENT2. The solid and dashed lines represent the predictions for t^ v with uncertainties based on 
the results in Eqs. dij) ) which are obtained from t2, y for y = r, n,pi (dots with error bars at In a = 0). All results are for 
n f = 5. 



The results for t2, Ta and ti. Pa with a = 2, 3, 5, 10 and 1/2, 1/3, 1/5, 1/10 for nf — 5 arc displayed in Fig. @, where 
we use the type of presentation from the last column in Fig. |l|. Comparing to the results displayed in Fig. |j] we see 
that for increasing values of In 2 (a) the asymptotic regime is shifted towards smaller values of f2. Thus for estimating 
the values for i2,r a and t2 tPa we now use the average of the lowest five bins with log 10 Q < —3.75 and adopt the error 
at log 10 Q = —3.75 as the uncertainty. The results are displayed as horizontal lines in Fig. [| and also summarized in 
Fig. |^. In Fig. |we have also shown the theoretical predictions for £2,t q and t%, Pa based on Eqs. ( |37| ) and the values 
of S\ and s 2 from Eq.(^5|) with iif = 5 for their respective central values (solid lines) and for their one-standard 
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5 (using a presentation analogous to the 



r,r a ,p a , for n; 

one in the right column of Fig. hi) separated according the color factor contributions proportional to CaCp (blue), Cp (red) 
and CfTu; (green). In the limit tt — » the Cp contributions are zero due to nonabelian exponentiation, and no error estimate 
is given. 



deviations (dashed lines) . The agreement is excellent and reassures the form of Eq. (|3l]) . 

Contributions from Different Color Factors. EVENT2 can give the 0{a 2 s ) contributions for the distributions separated 
with respect to the color factors C f ,CaCf and CpTnf and it is thus straightforward to determine the color factor 
components of the constants si and S2- The numerical determination of the Cp contributions is particular important 
since due to the nonabelian exponentiation property of the soft function si and S2 do not have any term proportional 
to Cp. 

Figure |B| shows the results from last column in Fig. ^ for i 2 . r , t2. Ta and t 2tPa separated according to the three types 
of color factors. 4 It is conspicuous that the different color factor components approach their asymptotic value at 
il = at much smaller values for fi as compared to the sum of all color factors shown in Fig. ^. In particular, for the 
Cp contributions the compatibility with zero becomes only apparent for log 10 f2 « —4. Thus our results for the color 
factor components have somewhat larger errors than for their sum that was obtained above. To obtain our results we 
use the average of the lowest five bins with log 10 Q < —3.75 and adopt the error at log 10 f2 = —3.75 as the uncertainty. 
For the color factor contributions of s± and s 2 we find 



S X,C% = S 2,C% = , 

sx,c A c F = -117.6 ±4.5, 

73.0 ±5.1 (ti/ = 5) 
59.5 ±5.1 (n f =4) 



S\ ,71 f 



S2,C A C F 
S2,n f = 



= -12.3 ±2.6, 

-2.2 ±2.8 (n f = 5) 
-1.8 ±2.8 (n f = 4) 



(46) 



where we have adopted zero for the Cp contributions as required by the nonabelian exponentiation property. As for 
our results for the sum of all color factors, we have checked that the results in Eq. (Mm arc fully compatible with 



the color factor contributions of i2,r Q and t2 lPa obtained from EVENT2 for a = 2, 3, 5, 10 and 1/2, 1/3, 1/5, 1/10. A 



determination of the color factor components of s\ for was also carried out in Rcf. |26 . Transferred to our notation 



4 Since we found numerical instabilities for the CpTrif contributions obtained from EVENT2 for log 10 f2 <^ —3.25 we determined these 
contributions from subtracting the Cp and CaCp terms from the results for the sum of all color factors. 
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their result reads s 1C 2 = 8.3 ± 1.8, Sx,c A C P = —120.0 ± 2.0, si,„ /= 5 = 71.7 ± 1.7 and si,„ /=4 = 57.3 ± 1.3. For 
the Cp term a quite small error is claimed, rendering the result incompatible with zero. This is because in Ref. ]2(| 
EVENT2 was run with linear binning, and only results for log 10 fl > —3 were used for the analysis assuming that the 
asymptotic values can be extrapolated from them. While our analysis validates this approach for the sum of all color 
factors, it can be seen from Fig. || that it fails for the Cp color factor contribution of si. For the CaCf and the n/ 
contributions of s\ the results of Rcf. |^6| are compatible with ours, but we believe that our error estimate is more 
appropriate. 



V. SOFT FUNCTION GAP 

The curves in Fig. |^ show a rather poor behavior of the shape of the soft function in perturbation theory. This 
behavior can be improved considerably when subtractions are applied to the partonic soft function S par t that remove 
the O(Aqcd) renormalon in the partonic threshold at = p"l|]. 5 This renormalon is very similar in nature (but 
independent of) the well-known O(Aqcd) pole mass renormalon in the heavy quark mass threshold in the partonic 
on-shell limit q 2 — m 2 = |2?], ^8|. To remove this renormalon order-by-order in the perturbative expansion it 
was suggested in Ref. |ll| to introduce the scale-independent gap parameter A in the soft function model function 
as shown in Eq. (|l2|). At this level A represents an additional model parameter that can compensate numerically 
for the divergent higher order behavior of the soft function caused by the O(Aqcd) renormalon. To obtain a gap 
parameter that is more stable in perturbation theory and allows a meaningful determination from experimental data 
it is mandatory to remove the renormalon contributions in the partonic soft function by explicit subtractions. This 
can be achieved by writing 



A = A + 8 A , 



(47) 



where <5A is a perturbative series 



SA = SAi + 8A 2 + 6A 3 + ... 



(48) 



that contains exactly the same O(Aqcd) renormalon as the soft function. Starting from Eq. ( |To| ) and shifting the 
integration variables the soft function can then be rewritten as 



dV 

— oo 



d£~ s p ^(i + ~e + -SA,e-~e--sA^) / exp (£+-A,r-A) . 



(49) 



To cancel the renormalon between the partonic soft function and the series SA order-by-order in the a s expansion 
we now have to expand Eq. (f49| ) in the <5A^ simultaneously with the expansion for the partonic soft function S'part = 

•Spart + ^part + ^part + • ■ •> S0 that 



"Spart 



SA,?) = 5° art (i ± ,/x) + 



?part(^,M) -<5A 1 (^ + -| r )5p a 1 ,(^, M ) 



"S'part 



d 

U+ ' dt 

^2 



+-^){sA 2 sl 1 . t (e ± ,ri + 5A 1 s 1 pait (e ± ^)} 



\de+ 2 + d£- 2+ de+de-J 2 ,Mj 



(50) 



where 5Ai and are of 0(a l s ). We stress that it is mandatory to use the same renormalization scale fi for the 

expansion of 5 A and S'part to achieve the renormalon cancellation, see e.g. Ref. ^9|. In Rcf. jll| a definition of the 
series in SA was proposed based on a ratio of moments of the partonic soft function with a finite cutoff. Since the 
soft function has an anomalous dimension, SA and therefore also A are /i-dependent. The moment definition does, 
however, not allow to formulate a consistent RG running of the gap parameter A due to logarithmic terms of arbitrary 



5 This renormalon in the soft function is the origin of the O(Aqcd) renormalon identified in Ref. |13| in the perturbative expansion of 
the thrust distribution in full QCD. 
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high powers in the evolution equations at high orders and because the resulting evolution equation is not transitive [p0[ 



Position space gap parameter. We can define a gap parameter with a consistent RG evolution from the position space 
soft function 



S(x u x 2 ,,u) = d£+d£-e- u+Xl e- u X2 S (£+,£-, n) 



d£ + d£~ e~ u Xl e~ a X2 



/+oo r + oo 

d£'+ / d£'~ S part {£+ - 5AJ- -f~ - SA, n) U P {£'+ - A,£'~ - A) 
-oc J — oo 



S par t(xi,x 2 , l i) \f CKp (x l7 x 2 )e- lA ^ +x ^ 



where 



and 



Axp(xi,a;a) - I d£ + d£~ e - a+x > e~ a ~ X2 f cxp (£ + ,£~) 



S paIt (xi,x 2 ,fi) = S pait (x 1 ,x 2 ,fi)e~ lSA( - x ^ . 
Since the function S par t(xi, x 2 , /x) is supposed to be free of the O(Aqcd) renormalon we can use the condition 



Re lE 



dh\{ix\) 

to derive an explicit expression for 6 A, 



lnS part (xi,x 2 ,LJ,) 



x 1 =x 2 = {iRe~'E)-i 



(51) 
(52) 
(53) 
(54) 



5A(R,fi) = Re lE 



\TiS paYt {x u x 2 ,ii) 



c?ln(i:ri) L r ~ 'J x 1 =x 2 =(iRe~iE)-i 
6&i(R,n) + SA 2 {R,fi) + 6A 3 (R,fi) 



(55) 



where the scale R is free parameter. This position space method was used before in Rcf. [|3Cf1 to derive a short-distance 
jet mass definition from the jet functions B± that appear in the factorization theorem (|6|) for massive quarks in the 
resonance region. The explicit results for the SA^s up to O(a^) read [L^r = ln(4j), a s = a s (/x)] 



<5A 2 (i?, M ) 
= Utt 



Utt 

Cts 

in 



7s ° + 2r°i Mii 



Utt 



2ti A) + 7] + 2 ( /3o 7 o + rj ) L^ R + 2(3 T° S fl 



808 22 



C A C F [- — --^+28C 3 +[- — + -^)L, R --L 



536 8 



9 



J HR 



224 8 2 160 32 2 

C F Tn/ ( — + -7T + — + y L^ R 



5A 3 (R,n) = (^) 3 [ 4 Sl p + 2t li 3 1 + 7s 2 + 2(4^/32+^7° + 2A)7 i 1 +r2)L /1 * 



Utt 



C A C F Tn f ( ^tt 



2(2/3 2 7s ° + ^1* + 2/3 riji2 i? + _^ r o L 3 fl 
62012 104 



81 ' 27 - 2 -^ 4 + 352 C3 )i^ 



14240 176 



27 



7T I LpR - 



3872 3 



32816 128 



it 2 \L 



^R 



C F Tn f [ 4tt 



C F (Tn,) 5 



/440 



V 3 

3200 128 



81 9 
128C 3 Ji / ut + 32l£ R 
1280 



9248 64 



27 9 



n 7 \L 



->liR 



2816 3 



81 27 



7T l-kfiii T^L,, n — 



27 



HR 



512 



( ^-C A - \Tn s 



(56) 
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The 0(al) corrections is only partially known as it depends on the unknown 3-loop non-cusp anomalous dimension 
7^ of the soft function. Using that A = A + 8 A is scale-independent it is straightforward to derive the anomalous 
dimension of the gap parameter A(R,fx), which turns out to be just proportional to the cusp anomalous dimension 
of the soft function in Eq. (jig), 



-—A(R,fi) = -——6&(R,n) = -2Rei*r.[a B (n)]. 
a m fj, a In /i 



(57) 



Up to 0(a^) this intriguing fact can derived explicitly from the results in Eqs. (56). To all orders it can be proven 
using the expression (22) for the position space soft function and the all-orders definition for K in Eq. (0). The 
solution of the RG equation reads 



A(R,fx) = A(R,n a )~Re^ E ij(T S: n,no 



(58) 



where u) is given in Eq. dT 



R-Evolution The definition of the gap in Eq. (|5q ) depends on the choice of the scale parameter R. Since SA(R, /x) 
represents an infrared subtraction of low-energy fluctuations from the partonic soft function, one can interpret the 
scale R as an infrared cutoff governing infrared fluctuations that are absorbed into the rcnormalon-frec gap parameter 
A(R, (j,). In this respect R differs significantly from the renormalization scale fi, which governs ultraviolet fluctuations 
that are absorbed into coupling constants. Like the renormalization scale /x, the scale R should be taken in the 
perturbative regime, to allow for a perturbative description of the evolution of A(R, /i) in R, but it should also be 
close to the typical scales that are governing the soft function. Moreover to avoid large logs [ijR should be of order 
one. For the description of the peak region for massless jet event shapes or the massive jet invariant mass resonance 
region one is sensitive to details of the shape of the soft function. So the typical range of scales for fi and R in this 
case is around 1 — 1.5 GeV. For describing these distributions in the tail away from the peak and resonance region 
one can use an operator product expansion and only global properties of the soft function in the form of moments are 
relevant. So here the perturbative contributions are dominated by larger scales, and for /i and R larger scales should 
be adopted as well. 6 

Besides these constrains the choices of fi and R are arbitrary. It is therefore useful to also have an evolution equation 
for the gap parameter with respect to R. A detailed study of the required formalism for evolving A(i?, R) in R and 
its relation to rcnormalons was carried out recently in Rcf. |31 and we refer the interested reader to this work. For 



the renormalon-frcc soft function gap parameter the devolution equation has the form 



dR 



A{R, R) 



R J2 



a s {R) 

47T 



R 



OsiR) 
47T 



To* 



OsiR) 
47T 



7? 



a s (R) 

47T 



if 



(59) 



where the first three terms can be obtained from Eq 
0, 

808 22 



7o* 



7i 



7* 



Ca Cf 



C\ Cf 



27 



9 



280 



/ 35552 662 



27 



1232 



3) and read 



C F TnJ™ 



V 81 



Cf {Tn } f 



'3584 128 



27 



V 81 



V 27 9 



C 3 )+C A C F Tn f 



AC F Tn f it 2 



■4si 



22784 524 



27 



81 



11 



-C A 



-Tn, 



448 



(60) 



At N fe LL order the analytic solution reads [tn = — 2ir / (f3oa s (R))] 



[A(R,R)-A(Ro,Ro)] 



N LL 



A 



QCD ^3 ( 1, 

3=0 



r(-&i-j,ift) - r(-6i-?,^ 



(61) 



Because the one-loop subtraction <5Ai term only contains a logarithm of fJ./R it is also required in practice to choose R < fi to have a 
subtraction with the correct sign at the one-loop level. 
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where [7« ee 7„V(2/3 ) 



50 = To , 

51 = 7 f -f&i+^To* 



5o 



7 2 * - f&i+MT* + 



l + bi)b 2 + \[bi + b. 



To* 



The terms 6j are the coefficient of the large-i expansion of the function b(t) = 1 + b±/t + b^jt 2 
the relation 



In — 

Ro 



«.(«) 



■s 



where the first few terms read 



A 

2/3 2 



6 1 - 7T^2> 



(iio) 

/?i 2 - P0P2 



dtb{t) = G{t Ro )-G(t R ). 



b, = 



4/3 4 



63 = 



/3? - Wft + /? 2 & 



and Aq^, d is the N fc LL order approximation of 



A 



QCD 



Re G{tn) = R a e G{tR o ] 



(62) 

that appears in 
(63) 

(64) 
(65) 



which is the familiar definition of Aqcd from the strong coupling. 

From the solutions for \i- and /devolution given in Eqs. ( [5cf ) and (61) one can relate A for arbitrary R and /i values 
through the relation 



A(R,fi) = A(Ro,v )-Re^ Q(T s ,Ro,fx ) +u>{T a ,n,R) 



Aqcd ^ Sj [T(-bi-j,t R ) -T^-bi-jM 

3=0 



(66) 



VI. NUMERICAL RESULTS 



A detailed analysis of the impact of the renormalon subtraction for the soft function has been given recently in 
Rcf. jllj using an approximation for the at that time unknown 0{p? s ) corrections. Since for the full 0(a 2 s ) soft function 
and the gap parameter obtained in this work we find qualitatively similar results, we do not repeat such an extended 
analysis here. In the following brief numerical analysis we intend to illustrate the impact for a few examples and to 
demonstrate the importance of accounting for the evolution of the gap parameter A(i?, /x) in (i?, /i)-space, which was 
not available for the gap parameter used in Ref. [ jTl| . 

In Fig. U the soft function S(£,£,fJ.) is plotted with and without renormalon subtraction for [i = 1.5 GcV and 
1.7 GeV (a s (1.5 GeV) = 0.3285) and different choices of R for the renormalon subtracted soft function. For all curves 
we have used nf — 5 and the soft function model <!5 m od with the parameters A = 0.55 GeV and (a, b) = (3.0, —0.5). 
For the soft function without renormalon subtraction (i.e. 5A = and A = A) we have used A = 0.1 GcV, and 
the tree-level (long-dashed black lines), 0(a s ) (dotted red lines) and 0(a 2 ) (dotted-dashed blue lines) results are 
displayed. The results show the rather poor perturbative behavior already observed before in Fig. ^| In particular, 
the soft functions at 0{a s ) and 0{a 2 ) have unphysical negative values for small values of I. For the soft function 
with renormalon subtraction we have used A(1.0 GeV, 1.5 GeV) = 0.1 GeV and computed A for the (R, /x) values in 
the different panels using Eq. (p6J) at NNLL order. 7 In the different panels the tree-level (solid black lines), 0(a s ) 
(lighter solid red lines) and 0(oq) (dashed blue lines) results are displayed. The results for the soft function with 



7 We have set the unknown non-cusp term 7^ to zero for the NNLL order R-evolution of A(_R, fj,). 
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FIG. 6: Dependence of the soft function S{i, t, /i) on the renormalization scale fi and the infrared subtraction scale R at 
tree-level, one-loop and two-loop with and without renormalon subtraction. In the middle column the solid and long-dashed 
black lines coincide. See the text for details. 
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FIG. 7: Location of the maximum of the soft function S(£, £, /J,) with renormalon subtraction as a function of R for /j, = 1.5 and 
1.7 GeV. The solid curves account for the correct evolution of the gap parameter A(R,/j,), while for the dashed curves a fixed 
gap parameter A = 0.1 GeV is used for all values of R and fj,. 



renormalon subtractions show a substantially improved perturbative behavior. Concerning the shape and the peak 
location, the 0{a 2 s ) corrections are particularly small. It is also conspicuous that the curves do not have unphysical 
negative values except for (R,n) = (1.5, 1.2) GeV and (1.7, 1.4) GeV where the difference of R and fj, is small so that 
the subtractions at 0(a s ) are not sufficiently large (see footnote 1). 

For the stabilization of the peak position with respect to variations of R it is important to use the consistent evolution 
of A in the (R, /i)-plane. In Fig. f^the position of the maximum of the C(a|) soft function with renormalon subtractions 
is displayed as a function of R for the soft function model used also in Fig. y and for fi = 1.5 GeV and fj, = 1.7 GeV. The 
solid lines show the peak location with a consistent choice of A(i?, fj) using A(1.0 GeV, 1.5 GeV) = 0.1 and Eq. ( |66| ) 
to obtain A for other (R, /i)-values. The dashed lines, on the other hand, show the results when A = 0.1 GeV for 
any R and fi. The curves show a considerably weaker i?-dependence when the correct evolution of the gap parameter 
A(R, /j,) is accounted for. Although this is expected due to the fact that A = A + <5A is i?-independent, it is reassuring 
that this is also reflected in the final result for the soft function given that i5A is implemented by subtractions that 
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change the shape of the soft function substantially, whereas A is implemented by a simple shift of variable £. Note 
that A = A + (5A is also /t- independent. However, since the soft function has a non-zero anomalous dimension, it 
is not expected that accounting for the proper //-evolution of A also leads to a smaller //-dependence of the peak 
location. This is confirmed by the vertical separation of the two solid curves that is in general not smaller than for the 
corresponding two dashed curves. However, it is clearly visible that the vertical separation for the solid curves (which 
account for the consistent evolution of A) is approximately i?-independent as compared to the vertical separation 
of the dashed curves (which have a fixed value of A). This shows that the variation of the location where the soft 
function has its maximum with changes of \i is only i?-independent when the proper evolution of A is included. 
This emphasizes once more the importance of having a consistent formulation of the R- and /x-evolutions of the gap 
parameter A(R, /i). 



VII. CONCLUSIONS 

In this work we have determined the partonic hemisphere soft function S , pa rt(^ + , , lA at 0(a 2 ). The partonic soft 
function S vax t is an essential ingredient for the theoretical description at NNLL order of thrust and heavy jet mass 
distributions in the dijet and resonance limit, where perturbative and nonperturbative contribution can be separated 
according to the factorization theorems (||) and (^). Using general properties of the partonic soft function originating 
from its renormalization group structure and its exponentiation properties, we have been able to determine the 0(a 2 s ) 
corrections up to two constants which were determined numerically using information from the MC program EVENT2 
by Catani and Seymour. 

To remove the O(Aqcd) renormalon of the partonic threshold in 5p ar t at £ = one can implement a subtraction 
procedure devised in Ref. |Tl] ] that is based on a gap parameter. Using this gap subtraction scheme one can remove 
order-by-ordcr this renormalon contribution in the perturbative parts of the factorization theorem that otherwise 
causes instabilities for numerical determinations of the first power correction from numerical data. From the result for 
the partonic soft function in position space representation we have defined such a gap subtraction scheme which has 
the virtue of having a consistent renormalization group evolution in the renormalization scale /t and the subtraction 
scale R. This property is important for having a coherent theoretical description of event shape distributions in the 
peak and the tail region. 
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APPENDIX A: HARD COEFFICIENT, JET FUNCTION AND CONSISTENCY CONDITION 



Hard Coefficient Hq. The 0(a 2 s ) fixed-order expression of the hard coefficient in the dijet factorization theorem (||) 
for massless jets can be derived from the results in Ref. |l6|, |l7|, |l8| and has the form [a s = a s (/t), Lq = ln(Q 2 //i 2 )] 

C F a s 



Hq{Q,v) = 1 



5J l<s 

+ Ca Cf 
+ C F Tn f 
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22,3 
9 Q 
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81 27 9 s 
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The RG evolution is described by Hq(Q, /i) = Hq(Q, (1q)Uh q (Q, Ho, m)j where 



Uh q {Q,Ho,h) = CX P 



w(r c , n, no) In 

The functions <D and if are given in Eqs. ( |17| ) and we have 

it = -rLp. (i = o,i 5 ...), 



2 

Q 2 



2X(r c ,7 c ,^,/i ) 



(A2) 



7c ° = -6C/T, 

7, 1 = Cj ( - 3 + 4tt 2 - 48C 3 ) + C A C F U^. ^ 



52C; 



27 3 

where the terms in the cusp anomalous dimension T c are given in Eq. (fL9| 



260 4 

Gp i Tit 7T 

1 27 3 



(A3) 



Jet Function J . The 0(a 2 s ) fixed-order expression for the jet function for massless quarks was computed in Ref. |L9| |, 
where the result was given in Laplace space. In momentum space the jet function reads 

[a. = aM Cl + = l/^[0( S )ln"( s /^)/( s / M 2 )] + ] 



J(s, M ) = 5(a) + (^-) - 3 < + + ( 7 - O 



(A4) 
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and its renormalization group evolution reads J(s,/x) = J ds'Uj(s — s', jj,, [io)J( s 'i /j,q) with 



Uj(s,fi,fx ) 



Mo r (-^( r J>M!Mo)) 



2u + i<i(r J , Al , A t ) , 



(A5) 



where the functions to and K are defined in Eqs. (p"7|), the cusp and non-cusp anomalous dimensions are given in 
Eqs. dAS| ) and the definition of the plus function given in Eq. ([30j) is used. In position space the jet function has the 
form [L£ = ln(i£e 7E ^ 2 )] 
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and its renormalization group evolution reads J(£,/z) = £//(£, /x, /Zo) J(£, ^0) with 

1 



Uj(£, fJ,, no) = exp 



(T j, /x, /x ) In u£e 7£ > + K(Tj,jj,ii,ho) 
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where the functions Co and K are given in Eqs. (jl7j). The cusp and non-cusp anomalous dimensions read 

r^=2r* usp , (i = o,i,...), 
7°/ =6C F) 

„i_^/ Q , AQ ,\ , A769 22 , Qn ,\,^^/ 484 8 , 



7.) = C F (3 - 4tH + 48C 3 J + C A C F + -tH - 80C 3 j+C F Tn f {-— - -tH j , (A8) 

where the terms in the cusp anomalous dimension can be read of from Eq. (p"9|). 

Consistency Condition and Soft Function Evolution Factor U s . In Rcf. ^ it was shown that the renormalization 
group evolution [/-functions of the hard coefficient Hq, the jet function J and the soft function S in the factorization 
theorem (^|) for jets initiated by massless quarks are related by a consistency condition. It can be used to determine 
the renormalization group evolution functions of the soft function. In momentum space the consistency condition 
reads 



, fi, im>) = Q ^U Hq (Q, fi, ft ) U J (Qe ± , fi , fj) , (A9) 

and in position space it has the form 



U a (xi i2 , fi, Mo) = \J U Hq (Q, fi, fio) Uj (^pp fio, fi) ■ (A10) 

Using the relations 

ui(T + T', fi, fio) = Q(T, fi, fi ) + Q(T', fi, fi ) , (All) 

K(T + V, 7 + l', ^ fio) = K(T, 7, fi, fi ) + K(V, 7', ft, fio) , 
w(r,/x,/x') = - u>(T,fi',fi) , 

K(-T,^,fi ,fi) - (DQ,/io,m) = - K(-T,^f,fi,fi ) +Q(^,fi,fi ) - ij(T,fi, fio) In , 
it is then straightforward to derive Eqs. (O) and (J29|). 



Hard Coefficient H m . Recently the jet function B± was determined at 0{a 2 s ) in Ref. |3(J]. This jet function describes 
low-energy invariant mass fluctuations of jets initiated by massive quarks in the factorization theorem ([|). Using 
the result for the NNLL anomalous dimension for B± and the soft function S, and the consistency condition for the 
/i-evolution of the contributions in the factorization theorem for massive jets in Eq. (|^), one can derive the NNLL 
anomalous dimension of the hard coefficient H m in Eq. (^|). For the results for the 0{a 2 s ) jet function B± in momentum 
and position space, which match exactly the conventions used here, we refer to Ref. fl30[ . The computation for the 
evolution factor for H m is most conveniently done in position space representation. We start from the relation between 
the position and momentum space jet function, 



B±(x,fi) = ldse- iys B±(s,fi). (A12) 
The RG evolution of the position space jet function is described by B±(y, fi) = Ub(u, fi, fio)B±{y, fio), where 



UB(y,fi,fio) = exp 



u)(T B ,fj,,fi )hi (iye lE ^ + K(T b ,Jb, fi, fio) 



(A13) 



where the functions Co and K are given in Eqs. (|l7|). The cusp and non-cusp anomalous dimensions read 

T'b = Tcusp j (i = 0, 1, . . .) , 

"fi =ic F , 

,h = A C F _ |„s _ 20<J ) + Cf Tn , (i^ - f ) . (AH) 
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The renormalization group evolution of the hard coefficient H m is described by H m (m,,Q/m, fj, m , /i) = 
H m (m, Q/m, /io)[^F/ m (Q/m, Mo> m) and using the consistency condition the evolution factor Un m can be related to 
Ub and U s , 
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Q 



\m ) L Vm /J L 

Using the results for Ub and U s in Eqs. (A13) and ( |T6| ) and the relations (All) we obtain 
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where the cusp and non-cusp anomalous dimensions read 
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APPENDIX B: FIXED-ORDER r Q AND p a DISTRIBUTIONS 

At 0{a 2 s ) the fixed-order cumulative r Q distribution for fi = Q reads 
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At 0(a 2 a ) the fixed-order cumulative p a distribution for fx = Q reads 
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From the cumulative distributions in Eqs. (Bl) and (B3) one obtains the unintegrated distributions (\ / cr^)da d ^ et / dr a 
and (1 / <Jo)d(T dlsct / dp a by differentiation. This entails applying the following replacement rules [y — r a or p a , k = ^t 2 ]: 
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APPENDIX C: FOURIER TRANSFORM 



Given a momentum space variable t with mass-dimension j the Fourier transform of delta-functions and plus- 
distributions of t into position space with the variable y can be derived from the relations 

dte ~' tvm \ujy =r(e) • (C1) 



and 



n=0 
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From this result one finds [£™ + = l/n j [6(t) hi n (t/^)/(t/fi j )] + , L y = In^ye™ /j, j )}] 

fdte-^5(t) =1, [die-** = -^-y^-^Ca, 

Jdte-*»# t+ = l -L y + ^. 
The transformation from position back to momentum space can be derived from 
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which leads to 
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The cumulative distributions can be obtained directly from the position space distributions using the relation 
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This leads to [L yK = ln(i?/e 7E K), Lfi K = ln(0/n) 



4/ 

27r i(y — iO) 

dy 
— e 
2tt i(y — iO) 

e in y L l* 

2t: i(y — iO) 



6(0), 

-9(0)L Qk , 
6(0) 



A more general formula relevant for renormalization group improved computations reads 
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